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In this paper we present a Monte Carlo calculation of the critical temperature and other ther- 
modynamic quantities for the unitary Fermi gas with a population imbalance (unequal number of 
fermions in the two spin components). We describe an improved worm type algorithm that is less 
prone to autocorrelations than the previously available methods and show how this algorithm can 
be applied to simulate the unitary Fermi gas in presence of a small imbalance. Our data indicates 
that the critical temperature remains almost constant for small imbalances h — A/j,/ef ~ 0.2. We 
obtain the continuum result T c = 0.171(5)£f in units of Fermi energy and derive a lower bound on 
the deviation of the critical temperature from the balanced limit, T c (h) — T c (0) > —0.5eFh 2 . Using 
an additional assumption a tighter lower bound can be obtained. We also calculate the energy per 
particle and the chemical potential in the balanced and imbalanced cases. 



I. INTRODUCTION 

The Fermi gas at unitarity - a dilute system of two- 
component fermions interacting with divergent scattering 
length - is a particularly interesting example of a strongly 
interacting fermionic system jl] |3| . In this case the den- 
sity sets the only relevant length scale and the system 
exhibits universal behaviour. At a certain critical tem- 
perature a phase transition into a superfluid state takes 
place. In contrast to the weakly interacting BCS limit 
the critical temperature at unitarity is of order of the 
Fermi temperature, and hence accessible for experimental 
study. Examining the mechanism behind this phase tran- 
sition promises valuable insights into high-temperature 
superfluidity. 

In the strongly interacting limit, perturbative methods 
are inapplicable and mean-field approaches involve un- 
controlled approximations, since the contribution from 
fluctuations becomes significant Hence, numerical 
methods have found wide application. The method un- 
derlying our study is the Diagrammatic Determinant 
Monte Carlo (DDMC) algorithm j|, which was devel- 
oped for the calculation of the critical temperature and 
was specifically designed to take advantage of the physi- 
cal properties in the unitarity limit. Here we will describe 
our tests and implementation of this algorithm and sug- 
gest several modifications that significantly reduce auto- 
correlations and therefore increase the efficiency. Prelim- 
inary results were presented in || . 

So far, most numerical studies were limited to the bal- 
anced case, when the number of fermions in the two 
spin components is equal. An imbalance in the particle 
number results in new interesting effects, which makes a 
detailed numerical study desirable ||, [7]]. On the other 
hand, an imbalance also leads to algorithmic difficulties, 
which are due to fermionic statistics. The Feynman di- 
agrams in the expansion of the partition function have 
different signs, depending on the number of fermionic 
loops. In the balanced case this problem can be avoided, 
since the diagrams of each order can be represented as a 
square of a matrix determinant B . In the presence of an 



imbalance this is no longer the case, so that the partition 
function cannot be used as a probability distribution for 
Monte Carlo sampling. A similar problem arises in lattice 
QCD, where a non-zero chemical potential renders the 
fermionic determinant complex. Many techniques deal- 
ing with the sign problem have been developed, the most 
straightforward of which is the "sign quenched method" 
Q, which we will make use of in our study. We found 
that at unitarity the sign problem is mild, which is also 
consistent with the observation that the superfluid state 
remains remarkably stable in response to increasing im- 
balance. 

In this work we present a numerical calculation of the 
critical temperature and other thermodynamic observ- 
ables for several values of imbalance. We begin with a 
review of the Fermi-Hubbard model an d the finite tem- 
perature formalism in Sec. Q. In Sec. Ill we define the or- 
der parameter for the phase transition and show how the 
critical temperature can be extracted from the numerical 
data. We then summarise our version of the worm algo- 
rithm in Sec. IV and explain how it can be generalised 
to the imbalanced case. Finally, our results for the bal- 
anced as well as the imbalanced gases are presented and 
discussed in Sec. [v[ 



II. FERMI-HUBBARD MODEL AT FINITE 
TEMPERATURE 

The Fermi-Hubbard model is the simplest lattice 
model for two-particle scattering. Its Hamiltonian in the 
grand canonical ensemble is given by 
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where £k = ;^Xw=i(l — cosfcj) is the discrete disper- 
sion relation, and cj^ (ck<r) the time-dependent fermionic 
creation (annihilation) operator. We set % = kg = 1 
throughout. We have chosen this simple dispersion re- 
lation for a better comparison with reference [0 where 
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the same relation was used. It is possible to speed the 
approach to the continuum limit by choosing a more com- 
plex dispersion relation . This will be explored in fu- 
ture work. This model describes non-relativistic fermions 
of two species labelled by a (which we will call "spin up" 
and "spin down") with equal particle mass m. The at- 
tractive contact interaction is characterised by the cou- 
pling constant U < 0. This coupling can be tuned so that 
the scattering length takes infinite value, by solving the 
two-body problem in the same way as it was done in [|] . 
The corresponding value is U = —7.914, in units where 
m = 1/2. We work on a 3D simple cubic spatial lattice 
with L 3 sites, periodic boundary conditions and lattice 
spacing set to unity. The continuum limit of this model 
can be taken by extrapolation to vanishing filling factor 
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To study the finite temperature behaviour we work 
with the grand canonical partition function in the imag- 
inary time interaction picture, Z = Tre - ^^, where (3 is 
the inverse temperature. The imaginary time direction 
remains continuous. Using Dyson's formula and expand- 
ing Z in powers of Hi generates a series of Feynman 
diagrams, where each 4-point vertex has one incoming 
line of each spin and one outgoing line of each spin. The 
Feynman rules assign a factor of (— U) to a vertex and a 
line represents a free (finite temperature) single-particle 
propagator, 
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where T T denotes the imaginary time ordering operator. 
The explicit form of the propagator in momentum space 
is given by @ 
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where = (1 + e^' £k Mct ') 1 is the occupation of 

the state (k, a) for free fermions. Additionally, each 
fermionic loop contributes a minus sign, with the con- 
sequence that the diagrams in the series have different 
signs. Since we are ultimately interested in thermal ex- 
pectation values of operators and thermal averages are 
calculated using the expansion of the partition function, 
it would be convenient to use this expansion as a prob- 
ability distribution to generate configurations for Monte 
Carlo sampling. For this purpose we need to rewrite the 
series as a sum of positive terms only. It was shown in 
Q that the partition function can be written as 



Z = ^2(-U) p det A t (S , p )det A±(S P ), 



(5) 



where S p denotes a vertex configuration (the spacetime 
positions of all vertices) and the matrix entries are the 
propagators Af,(S p ) = G^ r 0) (x l - x^r, - Tj), given by 
Eqs. (^) and (jj). Note that within this formalism the 
only degrees of freedom are the vertex coordinates and 



hence it is not necessary to distinguish between the dif- 
ferent ways of connecting them. If the chemical potential 
is equal for spin up and spin down fermions (the balanced 
case) we have det A^ det A^ = | det A| 2 , so that all terms 
in the series are positive. 



III. ORDER PARAMETER AND FINITE-SIZE 
SCALING 

The physical observable in the focus of our study 
is an order parameter for the phase transition to su- 
perfluidity To define the order parameter, which is 
related to the density of the condensate, we first in- 
troduce the pair creation and annihilation operators 
Pt(x',r') = 4, t (r')4, 4 (r') and P(x,r) = c^t^t). 
At the critical point the correlation function 



G 2 (xr;xV) = (T T P(x,r)Pt( x ' )T ')) 



(6) 



= iTr[T T P(x,r)Pt( x ' )T ')e-^] (7) 



is proportional to |x — x 



H-(i+v) 



as |x — x'| — > oo (in three 



dimensions), where 77 « 0.038 [|l2| is the anomalous 
dimension for the U(l) universality class. Hence, if no 
corrections due to irrelevant operators were present, the 
rescalcd integrated correlation function 



R(L,T) = L 1+r >(j3L 3 )- 2 V / dr dT'G 2 (xT;xV 
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would be independent of lattice size L at the critical j3 c — 
l/T c §]. In this case, all R(L, T) curves for different 
values of L would cross in a single point. Taking scaling 
violations into account, the function P(L,T) near the 
critical point can be written as a product of a universal 
analytic scaling function f(x), where x — L 1 /^^ and a 
correction term due to finite lattice size, 



R(L,T) = /(i x /^t)(l + c L 



(9) 



Here t=(T — T c )/T Cl c is a non-universal constant, and 
the critical exponents w w 0.8 and v% ~ 0.67 can be de- 
termined to high precision with various methods, see e.g. 
Jl2], 13 . Near the critical point we can expand Eq. (g) 



and keeping only terms linear in t we obtain 

R(L, T) = (/o + /i(T - T C )L 1 '^ + . . .)(1 + cL-" + ...). 

(10) 



Previous work fi,[l4| used a two-step procedure for deter- 
mining T c from this equation. First the crossings of the 
i?(L, T) curves for each pair of lattice sizes were deter- 
mined individually. By equating R(Li,Tij) = R(Lj,Tij) 
and using Eq. (10), a relation between T c and the crossing 
temperatures Ty can be derived, 



Tij -T c = ng(L l: Lj), 



(11) 
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FIG. 1. The relative difference 

(g(Li,Lj) — g(Li, Lj)) /g(Li, Lj) as a function of c (in 
lattice units), for several values of Li and Lj ranging between 
10 and 16. 
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FIG. 2. (Color online) A typical fit of the rescaled correlation 
function R(L,T) according to Eq. (px|). Here data was taken 
at four different lattices sizes and temperatures. For this fit 
X 2 /d.o.f= 1.4. The value for c was found to be -1.4(5). All 
quantities are given in lattice units. 



where 
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and k = c/o//i is a non- universal constant. In the refer- 
ences [Q and [0 the second term in the denominator of 
g (L i} Lj) is neglected, so that the second step of the pro- 
cedure simplifies to a linear fit of the T5y to the function 
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The critical temperature is the intercept of this fit. This 
simplification is justified if the constant c is sufficiently 
small. But already if c assumes values of order of unity 
the systematic error associated with this approximation 
can reach up to 20%, as shown in Fig. [t]. To avoid this 
systematic uncertainty we choose a different procedure 
for extracting T c from the numerical data. In our analy- 
sis we use Eq. @ directly to fit all data triplets (R, L, T) 
to a single function. An example of such a non-linear fit 
is shown in Fig. |[ This procedure has several advan- 
tages. Firstly, the previously described systematic error 
is no longer present, which can become relevant, since 
we found |c| > 1 in several cases (see Sec. |V|). Secondly, 
all information obtained from the simulations is used for 
the data analysis. In the original two-step procedure the 
(R, T) tuples for each L were fitted to a line separately, 
which involved two unknown parameters for each value 
of L. After the crossings of these lines were determined, 
the information about the values of R could no longer be 
used for the next stage of the analysis. From the cross- 
ings another linear fit involving two unknown parameters 
had to be made. For our example from Fig. ^| with 16 
datapoints, the original procedure would require four in- 
dependent linear fits (8 parameters) and another linear 



fit into which the errors of the previous fits propagate. 
The new method suggested here only requires a single 
non-linear fit of 4 parameters: /o, /1, c and T c , of which 
only T c is of interest here. 

The value for T c is obtained in lattice units and needs 
to be translated into physical units. Since the only phys- 
ical length scale at unitarity is determined by the den- 
sity, the corresponding physical quantity has to be T c /sf, 
where the Fermi energy is defined as Ef — (37r 2 ^) 2 / 3 . In 
the grand canonical ensemble the chemical potential is 
fixed and the corresponding filling factor v is measured 
for different values of lattice size. For sufficiently large 
lattices the values v(L) scale linearly with 1/L and an 
extrapolation to 1/L — > will yield the thermodynamic 
limit for the filling factor at a given chemical potential. 
Finally, the continuum limit for the critical temperature 
is taken by extrapolating to v — > j|, |l5| . 



IV. IMPLEMENTING THE ALGORITHM 

A. Balanced case 

The configuration space of diagrams can be sampled 
via a Monte Carlo Markov chain process: in each step one 
of the possible updates to another vertex configuration 
is proposed with probability W(S P — > S') and accepted 
with probability P(S P — > S' q ) = min(l,7t), given by the 
detailed balance equation 

nw(s p -> s' q )v^(s p ) - w(s p <- s' q )v^(s' q ), (u) 

where V^ Z \S P ) = {-U) p \ dot A(5 P )| 2 stands for the dia- 
gram corresponding to the vertex configuration S p . The 
requirements of detailed balance and ergodicity ensure 
that the configurations produced are indeed distributed 
according to the correct thermal probability distribution 
P z(S P ) = ±(-Uy\detA(S p )\ 2 . 



4 



The Monte Carlo estimator for a generic thermody- 
namic observable (X) — ^Tr[Xe~P H ] can also be found 
easily. If we denote the diagrams in the expansion of 
TrlXe^ 1311 ) by T>^ X \S P ), as we did for the diagrams in 
the expansion of the partition function, we can write 



(X) 



z 4^ 



[X] 



(Sp) 



(15) 

(Q {x > z \s p )) pz . (16) 



is the desired Monte Carlo 



Here QM{S P ) = 

estimator, given by the ratio of the weights, and (. . ) Pz 
stands for averaging over a sequence of Monte Carlo ver- 
tex configurations S p created according to the probability 
distribution pz(S p ). 

The diagrammatic expansion of the correlation func- 
tion Tr[T r P(x, r)Pt(x', T f )e-P H ] is similar to that of the 
partition function Z , but contains an additional pair of 2- 
point vertices at (x, r) and (x', r'). It is thus of advantage 
to sample these two series in the same simulation. In ad- 
dition to sampling the regular 4-point diagrams we allow 
updates that insert the pair of 2-point vertices ("worm 
vertices" ) into the configuration space [|J . In matrix no- 
tation this means that if the worm vertices are present, 
each of the Green's function matrices A CT gets an addi- 
tional row and column, coming from contractions with P^ 
and P respectively. We denote the extended probability 
distribution by pw(S p ) — -^\— U) p \ det A(S P )\ 2 , where 
the space of possible vertex configurations S p has been 
enlarged. This yields a change in the normalisation con- 
stant, such that the (extended) partition function now 
takes the form 



Z w = Z 1 



dr dT'G 2 (xr; xV) \ 



(17) 



where £ is an arbitrary parameter. The advantage is 
that the Monte Carlo estimator for the order parameter 
now becomes very simple: it is just a constant times the 
ratio of configurations with and without worm vertices. 
Other physical observables like the number density or the 
energy are still only measured when the system is in the 
"physical sector", namely when the worm vertices are 
not present. Due to the extension of the domain and the 
corresponding change in the normalisation constant we 
need to introduce an additional rescaling: 



{ X ) = /q(x,z)\ = ^w/ Q (x,z w] 

\ I Pz Z \ 



(Q (x ' Zw) ) PW 



- (Q (z ' Zw) ) PW 

A detailed description of the individual updates can 
be found in the appendix of J|. The idea behind the 
worm algorithm is that at low densities the major contri- 
bution comes from multi-ladder diagrams, these are con- 
figurations where the vertices are arranged into several 
vertex chains. Proposing updates that favour the cre- 
ation of such vertex chains will lead to higher acceptance 



ratios and increase the efficiency of the simulation. At 
low densities the acceptance ratios of the worm updates 
are an order of magnitude higher than the acceptance ra- 
tios of the simple diagonal updates, in which the vertices 
are inserted or removed at random. We found however, 
that the worm type addition and removal updates from 
the original setup suffer from strong autocorrelations, so 
that even after many successful updates the configura- 
tion does not change significantly 5|. To illustrate this 
we compare the measurements of the interaction energy 
(which is proportional to the diagram order) in the worm 
setup and the diagonal setup in Fig. [|. Both simulations 
used the same parameters and a comparable number of 
MC steps. Figure |] shows the blocking analysis of the 
relative error for the same quantity. Blocking is a widely 
used technique to estimate the error of an autocorrelated 
measurement. The single data points are arranged con- 
secutively into blocks of equal size, and each block is 
replaced by the average of the measurements it contains. 
Then the error is calculated for the resulting blocked sys- 
tem in the usual way. If no autocorrelations are present, 
the error will be independent of the block size. In the 
presence of autocorrelations N consecutive data points 
fluctuate less than N independent measurements. Hence 
the error will increase with block size, until the block size 
reaches the autocorrelation length of the system. 

Because of the large errors due to autocorrelations the 
worm setup is effectively less efficient than the standard 
diagonal setup. For this reason in the present study 
we employ the conventional diagonal updates, together 
with the modified worm addition and removal updates, 
as proposed in ||. This setup combines the advantages 
of the diagonal setup (weak autocorrelations) with the 
ones of the worm setup (high acceptance ratios). Below 
is a summary of all updates used in our simulation. 
For the modified updates we also give the values of the 
acceptance ratios 1Z (the other acceptance ratios can 
be found in [|J). The correspondin g form ulae for the 



imbalanced case will be given in Sec. [IV B 



Updates only concerning the worm vertices: 

• Worm creation/annihilation: insert /remove 
the pair P(x, r), pt(x',T / ) into/from the config- 
uration. In our setup the distributions for P and 
Pt are independent: both are distributed uniformly 
over the lattice, so that W(S P -> S p ) = (/3L 3 )~ 2 , 
where S p stands for the configuration S p with the 
additional 2-point vertices. The authors of Q de- 
scribe a setup in which the vertex P is selected 
randomly, and the vertex pt is then chosen in a 
spacetime hypercube of given extent around P. To 
avoid autocorrelations that can be associated with 
this scheme we employ the independent setup. The 
acceptance ratio is then 



K 



det A(S P ) 



det A(S P ) 



(19) 
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FIG. 3. The first 100000 numerical measurements of the interaction energy (a measurement takes place every 100 MC steps) 
with the worm setup (left) and the diagonal setup (right). Strong autocorrelations are visible in the worm setup. 
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FIG. 4. The blocking error analysis of the interaction energy with the worm setup (left) and the diagonal setup (right). The 
blocked error is much higher in the worm setup and continues increasing even for large block sizes. 



• Worm shift: shift the P^(x',r') vertex to other 
coordinates. This update is equivalent to the worm 
shift update in Q and involves a shift to a nearest 
neighbour on the lattice and a time shift in some 
interval around the old coordinates. 

Updates of the regular 4-point vertices: 

adding/removing a 4-point vertex (changes the dia- 
gram order). 

• Diagonal version: add or remove a random ver- 
tex. This is the most basic setup for changing the 
diagram order, however at low densities the accep- 
tance ratios are very low. 

• Modified worm-type updates: 

— Choose a random 4-point vertex from the con- 
figuration (which will act as a worm for this 
step). 

— Addition: add another 4-point vertex on the 
same lattice site and in some time interval of 



length At around the worm. 

— Removal: remove the nearest neighbour of the 
worm vertex (implies that addition can only 
be accepted if the new vertex is the nearest 
neighbour of the worm). 

The probability density for the addition update is 
then W(S P — > S p+ i) = l/(pAr), where 1/p comes 
from selecting the worm and 1/Ar from choosing 
the new time coordinate. Analogously for the re- 
moval update W(S P <— S p+ i) = \/{jp + 1) and the 
acceptance ratio becomes 



n 



detA(5 p+ i) 



dct A(S P ) 



(-U)pAr 
p+1 



(20) 



The modified worm setup still prolongs existing vertex 
chains like the original worm setup, but autocorrelations 
are significantly reduced since the worm changes with 
every update. This new type of updates can only be em- 
ployed in addition to the regular diagonal addition and 
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FIG. 5. (Color online) The blocking error analysis of the in- 
teraction energy. Red circles correspond to the pure diagonal 
setup and blue squares to a combination of diagonal and mod- 
ified worm updates with equal probabilities for each kind of 
update. 



removal updates. It works regardless if the pair of 2- 
point vertices is present or not in the configuration (the 
original worm addition/removal updates can only take 
place when the 2-point vertices are present). The accep- 
tance rates for this update are comparable with those for 
the regular worm updates. To demonstrate the increase 
in efficiency we compare the diagonal and the modified 
worm setup at low density, when the acceptance rates 
of the diagonal updates are particularly poor. We again 
consider the blocking analysis of the relative error for 
the interaction energy. As Fig. ^| clearly shows, the au- 
tocorrelation length does not increase in presence of the 
modified worm updates. The blocked error in this case is 
significantly lower due to the increased acceptance rate. 

The structure of the updates requires the calculation 
of a matrix determinant of a large matrix (with rank p 
up to about 6000) in each MC step. Since only a few 
elements of the matrix change in the course of one up- 
date (at most one row and one column) a recalculation of 
the whole determinant from scratch is not necessary. In 
our implementation we make use of the fast matrix up- 
date formulae Q which decrease the number of required 
operations to order p 2 instead of order p 3 . 



B. Imbalanced case 

The original DDMC algorithm relies strongly on 
the assumption of equal densities of the two fermion 
species. This assumption allows us to write the par- 
tition function ^ as a sum of positive terms only, 
and consequently to use it as a probability distribu- 
tion for Monte Carlo sampling. To study the imbal- 
anced case m ^ fii a generalisation of the algorithm 
is necessary. Due to the sign problem the function 



p w (Sp) = -^-{-U) p det At(5 p )det A+(S P ) is no longer 
positive for all configurations S p and can thus not be used 
as a probability distribution. To deal with this problem 
we will make use of the "sign quenched method" , which 
is based on the "phase quenched method" known from 
lattice QCD ||. The idea is to write the function pw as 
a product of its modulus and its sign, 



Pw(S P ) = — (-U) p \detA^(S p )detA^S p )\sign(S p ), 

ZjW 



and to use the positive function 



(21) 



Pw(s P ) = -i-i-uy 



|detA t (5'p)detA^(S*p)| (22) 



as the new probability distribution. The factor Z' w en- 
sures normalisation. This reweighting implies another 
change in the Monte Carlo estimator for a generic ther- 
modynamic observable (X). When sampling according 
to the sign quenched probability distribution p' w (S p ), re- 
lation (1181) becomes 



<*> = 



Q {X ' Z - ] ) p , (Q^(S p )sign(S p )) p , 
" (Q(z^)(S p )sign(S p )) p , w 



(23) 

The Monte Carlo estimators Q remain unchanged, apart 
from a multiplication with ±1 depending on the rela- 
tive sign of the two matrix determinants det AJ(S P ) and 
det A^(Sp). This representation of a thermal average in 
terms of the new probability distribution is mathemati- 
cally equivalent to the usual thermal average. However, 
numerical errors can become very large if the expecta- 
tion value of the sign in the denominator is close to zero, 
as it happens for the expectation value of the phase in 
QCD. For the unitary Fermi gas the sign remains very 
close to unity for small imbalances, as shown in Fig. |^, 
so that sign quenching is applicable for imbalances up to 
approximately 0.2e^. The restricting factor that keeps 
us from reaching large imbalances is not the sign, but 
rather the fact that even large values of Ap do not nec- 
essarily lead to large differences in the filling factors of 
the two components and hence the physical value Ap/ ef 
still remains small. This method works best close to the 
balanced limit and can provide a useful tool to examine 
the trend of the critical temperature for small deviations 
from it. 

The worm updates and acceptance ratios now gener- 
alise straightforwardly to the imbalanced case. In all for- 
mulae we merely need to replace the terms | det A| 2 by 
| det A^ det A^\. A slight drawback is that we now need 
to keep in memory two large matrices instead of one and 
update each of these matrices separately. Also the rel- 
ative error of the sign adds to the relative error of each 
observable. 
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FIG. 6. Schematic plot of the average sign near the critical 
point as a function of imbalance. The shaded area covers the 
range of values the sign can take at different values of lattice 
size and chemical potential. The lower boundary of this area 
is the "worst-case" curve of the sign, corresponding to lowest 
densities and largest lattice sizes used. 
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T (v = 0) = 0.1 73(6)E: 
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FIG. 7. (Color online) The critical temperature versus filling 
factor for different values of the chemical potential. The con- 
tinuum limit corresponds to v — >■ 0. The linear extrapolation 
(solid line) of the seven data points at lowest filling factors 
(filled circles) yields a continuum value of T c /ef = 0.173(6). 
The dashed line corresponds to a quadratic fit through all 
data points. 



A. Balanced results 

Before we include the imbalanced data we first present 
our analysis of the data at zero imbalance, for compari- 
son with previous results from Q. We performed simu- 
lations at eight different values of the chemical potential, 
corresponding to eight different filling factors. The lat- 
tice sizes varied between 4 3 for the highest filling factor 
and 26 3 for the lowest, so that the volume range in phys- 
ical units was approximately constant. As discussed in 
jl| , for sufficiently small v the critical temperature scales 
linearly with 2/ 1 ' 3 . This behaviour is seen for z/ 1 / 3 < 0.75. 
Our results and the continuum extrapolation are shown 
in Fig. ^. A line was fitted through the seven points with 
i/ 1 / 3 < 0.75, resulting in T c /e F = 0.173(6) - 0.16(1)^/ 3 . 
The goodness of fit is x 2 /d-0.f = 0.39. For comparison we 
also fit a quadratic through all eight data points, result- 
ing in a continuum value of T c jep = 0.188(15), which 
is in excellent agreement with the linear extrapolation. 
This confirms that sub-leading corrections proportional 
to v 2 ' 3 can indeed be neglected for sufficiently small v. 
In Fig. ^| we show the results for the fit parameters c, /o 
and fiT c , according to Eq. (|To|) . These parameters are 
smooth functions of the filling factor. This data shows 
that the non-universal constant c does indeed take values 
of order unity and thus cannot be neglected. 

Our final result for the critical temperature in physical 
units is thus T c /ep — 0.173(6). This value is signifi- 
cantly higher than the previous result from Q, where 
T c /sf = 0.152(7). We will also make a comparison with 
other results available from the literature. In jl6j the 
result of [Q was found to be in agreement with a con- 
tinuous space-time DDMC method. The authors of p"4j ] 



found an upper bound of T c /ef ~ 0.15(1). They used 
an auxiliary field Monte Carlo approach and extracted 
the critical temperature from the finite-size scaling of 
the condensate fraction, using the same procedure as 
in Q]. The difference between our results might be at- 
tributed to the approximation made through this fitting 
method. A less recent result by the same group JlT], |lj] 
is T c /ef = 0.23(2). Through extrapolating Monte Carlo 
results of low-density neutron matter, the authors of |l|] 
found a value of T c /ep — 0.189(12) at unitarity. Their 
value agrees with our result within errors. There are also 
results obtained with the Restricted Path Integral Monte 
Carlo method ^0|, T c /ep ~ 0.245, and an upper bound 
of T c /ef < 0.14 obtained with a hybrid Monte Carlo 
method plj . Results obtained with an e-expansion are 
also available [2^] . For comparison, the critical tempera- 
ture in the BEC limit is T BE c = 0.218e F - 

The authors of |2^] conjecture that the leading order 
change of the critical temperature is linear in kpr e , where 
kp = y/sF is the Fermi wavevector and r e is the effective 
range of the potential (r e = —0.3056 in units of lattice 
spacing for the Fermi-Hubbard model with a model 
independent coefficient. Our result for the linear slope is 
AT c /ef = —0.16(1), which is larger in magnitude than 
the value from Q. 

A similar continuum extrapolation can be performed 
for other thermodynamic observables, like the energy per 
particle and the chemical potential. The corresponding 
data together with the fits is presented in Fig. ^. The to- 
tal energy is composed of the kinetic energy -Ekin and the 
interaction energy E- mt = (Hi). An explicit expression 
for the former can be obtained from the position space 
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FIG. 8. (Color online) The non-universal constant c (left) and the first two coefficients in the expansion of the universal scaling 
function f(x) (right) in lattice units versus filling factor, see Eq. ( |fo| ) with /(0) = fo and /'(0) = fiT c . 
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FIG. 9. (Color online) The energy per particle (left) and the chemical potential (right) versus filling factor at the critical point. 
The linear fits were performed for data at filling factors v 1 ^ < 0.75. Filled circles indicate data included in the fit and empty 
circles stand for data excluded from the fit, see text for discussion. 



picture, 




= 2(i 3 (64c x - 64c x+j )), for any j, a. (26) 

The factor L 3 in the last line comes from summation over 
all lattice sites and the other numerical factors are due 
to summation over j and a. Since 2(<4.c x ) corresponds 
to the filling factor, the kinetic energy per particle can 
be written as 

E kin /L 3 v = 6(l-2(cic x+j )/v). (27) 

The Monte Carlo estimator for the interaction energy is 
given by Q^ Hl \S p ) = —(3~ 1 p, as shown in 



The results shown in Fig. ^| are obtained at T c , but the 
temperature dependence of the chemical potential and 
the energy per particle was found to be very weak. Also 
these quantities showed almost no dependence on lattice 
size L. For the energy per particle we obtain the contin- 
uum value E/Nef = 0.276(14). In units of the ground 
state energy of the free gas, Epc — (3/5)Nef, our result 
is E/Epc = 0.46(2). Since the expression for the kinetic 
energy ( p7| ) involves a difference of two quantities of com- 
parable size, large fluctuations can occur, especially at 
low filling factors. For this reason measurements of the 
energy per particle at lowest filling factor could only be 
performed on lattices with size L < 14, which is smaller 
than the lattice sizes used for the measurement of the 
critical temperature. We include this point in the plot in 
Fig. [| (left), but exclude it from the linear fit. The good- 
ness of fit is x 2 /d.o.f. = 2.1. Our result shows excellent 
agreement with the value E/Efg = 0.45(1) at T c quoted 
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in @. The value quoted in g is E/Ne F = 0.31(1), 
which in units of the free ground state energy roughly 
corresponds to E/E F q = 0.52(2). 

For the chemical potential at T c we obtain the contin- 
uum value n/e F = 0.429(9) with x 2 / d -°- f - = 2 -8- Our 
result differs from the value \ije F — 0.493(14) quoted 
in ||], but is consistent with the value [ije F = 0.43(1) 
quoted in Q. 

Since the chemical potential and the energy are ex- 
pected to stay almost constant at temperatures below 
T c we also make a comparison to values from the liter- 
ature obtained at zero temperature. In the zero tem- 
perature limit the quantities (J,/s F and E/E F q are equal 
and Monte Carlo estimates range between approximately 
0.40(1) and 0.44(1) [|f-§8). Our value for the chemi- 
cal potential falls within this range, the value for the 
total energy is slightly higher, which is consistent with 
the fact that the energy must increase at finite tempera- 
ture. These numerical estimates are consistent with ex- 
periment |9|-§l|. 

Finally, we make a comparison with recent experimen- 
tal studies of the homogeneous unitary Fermi gas. A 
direct measurement of the critical temperature and the 
chemical potential of the uniform gas has been presented 
in p^| . Their experimental value T c /e F = 0.157(15) 
agrees well with our result. However, the value of the 
chemical potential at the critical point /j/ej? = 0.49(2) 
differs from our value. Another experimental determi- 
nation of the critical temperature and thermodynamic 
functions, including the energy and the chemical poten- 
tial, is described in Their values T c /e F = 0.17(1) 
and /i/s F — 0.43(1) at T c show excellent agreement 
with our results. Their result for the energy per par- 
ticle E/Ne F — 0.34(2) at T c is higher than our value. In 
another experimental work p4| an estimate for the crit- 
ical temperature at zero imbalance is extrapolated from 
data at higher values of imbalance. 



B. Imbalanced Results 

Now we will present our results for the imbalanced case 
M4- 7^ A*t- Data was taken at 25 points, of which 23 lie 
within the regime of linear scaling, v x l 3 < 0.75. Out 
of these 23 points 7 are at zero imbalance, as discussed 
in the previous section. The two most common ways of 
quantifying imbalance are either through the chemical 
potential difference Afi/e F — — n±\/e F , or through 
the relative density difference /S.v/v = \v^ — v\\/{v\ + V\)- 
For the values of imbalance considered in our study 
these two quantities are proportional to each other, with 
Av/v = 0.122(2)A/Ve_F, as illustrated in Fig. Iiol The 



relative density difference shows no dependence on lattice 
size (the L-dependencies of v and Ai/ cancel each other 
out), but considerable dependence on the temperature. 
Also since Az/ is a small quantity, numerical fluctuations 
can become significant. Since the chemical potential dif- 
ference is less prone to numerical errors, we will use it 




0.25 



FIG. 10. (Color online) Relation between the chemical po- 
tential difference and the relative density difference at T c . 



from now on to quantify imbalance. 

The critical temperature T c /e F is now a function of fill- 
ing factor v and imbalance h — Afi/s F (in the following 
we assume that all quantities are in physical units and 
do not always write the e F factors explicitly). We are ul- 
timately interested in the continuum limit corresponding 
to v — and want to perform the corresponding extrap- 
olation. To achieve this all numerical data is fitted to a 
three dimensional surface, where the following assump- 
tions are made for the form of the fitted function: 

• At fixed imbalance the critical temperature is 
a linear function of f 1 / 3 , with slope a(h): 
T c {v,h = const) = T c (h) + a(h)v^ 3 . This is a 
generalisation of the relation valid in the balanced 
case. 

• T c (h) and a(h) viewed as functions of the imbal- 
ance h are analytic and can thus be Taylor ex- 
panded. 

• Due to symmetry in h all odd powers in the Taylor 
expansions of T c (h) and a(h) have to vanish. 

• T c (h) must be a non-increasing function of h. 
Hence the fitted function takes the form 



T c {v, h) — T c (h) + a(h)i> 1 / 3 . 



(28) 



If we expand T c (h) and a(h) to leading order in h the 
fitted function becomes 

T c (u, h)=T + T 2 h 2 + (a + a 2 h 2 )iy 1/3 . (29) 

This requires a linear fit of four parameters. The best fit 
yields T = 0.171(5), a = -0.154(9), T 2 = 0.4 ±0.9 and 
a 2 = -0.7 ± 1.9 with x 2 / d -°- f -= 0-43. Note that the T 2 
value corresponding to the minimal \ 2 is positive, which 
is forbidden by physical arguments. The \ 2 function is 
very flat along the T 2 direction, so that forcing T 2 — 
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v 1 ' 3 

FIG. 11. (Color online) Projection of the data onto the (z/ 1//3 - 
T c ) plane. Red circles denote the balanced data and blue 
triangles data at non-zero imbalance. The line corresponds to 
the constant fit (|33|). Dashed lines denote the error margins. 

results in x 2 /d.o.f.= 0.44. From the error on T 2 we derive 
the lower bound T 2 > —0.5. The best fit values for To 
and «o are i n excellent agreement with the ones obtained 
from the fit of the balanced data only. 

The error on the best fit value for a 2 is very large and 
the fit is consistent with a 2 — 0. Hence we also perform 
a fit to the function 

T c (is : h)=T Q +T 2 h 2 +a v 1 / 3 : (30) 

where T c (h) has again been expanded to quadratic order 
and the function a{h) has been replaced by a constant 
aQ. The best fit is 

T c (v, h) = 0.171(5) + 0.07(ll)/i 2 - 0.155(8)^ 1/3 , (31) 

with x 2 /d.o.f.= 0.41. This x 2_vaiu e is even lower than 
for the previous fit, which means that the data justifies 
dropping the a 2 term. The best fit result is still con- 
sistent with T 2 = and leads to a much tighter lower 
bound T 2 > —0.04. The other parameters To and ao 
agree with the results from the previous fit and the fit of 
the balanced data. 

Since our results indicate that T c remains almost un- 
changed in response to a weak imbalance, we also perform 
a fit to constant T c (h) and a(h), 

T c (v,h)=T Q + a is 1 / 3 . (32) 

This is the same function as the one used in the balanced 
case and corresponds to a straight line fitted through the 
projection of all data points onto the (z^ 1/,3 -T c ) plane, see 
Fig. [ll]. The best fit is 

T c (v, h) = 0.1720(45) - 0.156(8)^ 1/3 , (33) 

with x 2 /d.o.f.= 0.41. Again the result agrees with the 
previous fits. 



0.05 0.1 0.15 0.2 0.25 

Au/e F 

FIG. 12. The continuum limit of the critical temperature as 
a function of imbalance. The solid line is the value obtained 
from the constant fit (^), the shaded area corresponds to 
one standard deviation. The dashed line is the lower bound 
obtained from fit (^) and the dot-dashed line is the tighter 
lower bound obtained from fit (j3o|). 
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FIG. 13. (Color online) Three dimensional plot of the critical 
temperature versus filling factor and imbalance. The surface 
corresponds to the constant fit (^). 



We also performed fits using the jackknife method and 
several robust fits. All results were consistent with the 
minimal \ 2 n ^s. Table | provides an overview of the 
results obtained with the different fit methods. The val- 
ues for the parameters Tq and oq were obtained with 
high accuracy and are all in excellent agreement with 
each other, independently of the form of the fit function. 
Depending on the model assumptions two lower bounds 
could be derived for the leading order deviation of the 
critical temperature from its balanced value. Figure [l2| 
shows these two bounds compared with the value in the 
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TABLE I. Comparison of fit parameters obtained by different fit methods. 



T 



T2 lower bound 



on 



Sao 



a 2 



5ao 



xVd.o-f. 



balanced data 
fit to eq. (S) 
fit to eq. (3C) 
fit to eq. (32) 



0.173 
0.171 
0.171 
0.172 



0.006 
0.005 
0.005 
0.0045 



0.4 
0.07 



-0.5 
-0.04 



-0.16 
-0.154 
-0.155 
-0.156 



0.01 
0.009 
0.008 
0.008 



-0.7 



1.9 



0.39 
0.43 
0.41 
0.41 




FIG. 14. (Color online) Projection of the data for the energy 
per particle onto the (v 1 ^ 3 -E) plane. Red circles denote the 
balanced data and blue triangles data at non-zero imbalance. 
For comparison the fit at zero imbalance is shown. The points 
at non-zero imbalance tend to lie above the balanced fit line. 



FIG. 15. (Color online) Projection of the data for the aver- 
age chemical potential onto the (1/ 1 / 3 -//) plane. Red circles 
denote the balanced data and blue triangles data at non-zero 
imbalance. The solid line corresponds to the constant fit and 
the dashed lines indicate the error margins. 



balanced case. A three dimensional plot of the data to- 
gether with a constant surface fit is presented in Fig. [T|. 
Experimental determinations of T c as a function of h are 
an area of active research, and much larger values of h 
can be reached Q . 

A similar analysis was performed for the energy per 
particle and the average chemical potential ^l/ef = |a*T + 
fil\/2ep. Since with increasing imbalance interactions 
become suppressed, we expect the absolute value of the 
interaction energy to decrease. This in turn means an 
increase of the total energy, since the interaction energy 
is negative. As we did for the critical temperature we fit 
the energy in units of Efq to the function 

E(u, h) = E + E 2 h 2 + + af } h 2 y /3 (34) 

and obtain the best fit parameters Eq = 0.440(15), 
4 B) = -0.17(3), E 2 = 3.4 ± 2.2 and a 2 E) = -3.1 ± 4.5, 
with x 2 /d.o.f.= 2.8. These results are consistent with 
the balanced fit. The leading coefficient E 2 — 3.4 ± 2.2 
is no longer consistent with zero. We also perform a fit 
to the function 

E{v, h) = E + E 2 h 2 + 4 B V /3 (35) 

and obtain the best fit result 

E{v, h) = 0.444(13) + 1.9(3)/i 2 - 0.18(2)^ 1/3 , (36) 



which agrees with the previous result. The x 2 /d.o.f.= 
2.7. Figure |lj shows the numerical data. 

The average chemical potential is not expected to de- 
pend on the imbalance. Hence we fit our data to the 
function 

fi{u,h) = ( uo + 4 A ' ) ^ 1/3 ( 37 ) 

and obtain 

h) = 0.429(7) - 0.27(l)i/ 1/3 , (38) 

with x 2 /d.o.f.= 1.1. This is in very good agreement with 
our balanced result. A plot of the data and the fit is in 
Fig. |TJ. 

VI. SUMMARY AND CONCLUSION 

In this paper we presented a Monte Carlo calculation 
of the critical temperature and other thermodynamic ob- 
servables of the unitary Fermi gas with equal and unequal 
chemical potentials for the two spin components. For 
our study we developed a modified version of the worm 
algorithm, which is less susceptible to autocorrelations 
than the previously available methods. This algorithm 
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can be generalised to the imbalanced case using the sign 
quenched method. We first calculated the value of the 
critical temperature using only data at zero imbalance 
and found T c = 0.173(6)£f, which is significantly higher 
than the result previously obtained with the worm al- 
gorithm (ij. One possible explanation is the difference 
in finite size analysis methods as described in Sec. III. 



Our results for the energy and the chemical potential are 
E = 0.46(2)£ FG and fi = 0.429(9)e F . 

In the imbalanced case we extracted the dependence 
of the critical temperature on the imbalance h = A/i/e F 
close to the balanced limit. Our analysis is consistent 
with T c /ef = const, for the range of imbalances con- 
sidered {h ;$ 0.2). The value at h = extracted from 
a quadratic fit of the balanced and imbalanced data 
was found to be T c (h = 0) = 0.171(5)e F , which agrees 
with the value obtained from the balanced data only. 
We further derived a lower bound on the leading or- 
der term in the expansion of the critical temperature 



T c (h) - T c (0) > -0.5e F . With the additional assump- 
tion that the linear dependence of T c /sf on v 1 ^ remains 
unchanged in the presence of a small imbalance a tighter 
lower bound of T c (h) - T c (0) > -0.04e F could be ob- 
tained. We also analysed the behaviour of the energy 
and the chemical potential in the presence of an imbal- 
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